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An expression for the spectrum of quantization error in a discrete-time system 
whose input is a sinusoid plus white Gaussian noise is derived. This quantization 
spectrum consists of two components: a white-noise floor and spurious harmon- 
ics. The dithering effect of the input Gaussian noise on both components of the 
spectrum is considered. Quantitative results in a discrete Fourier transform (DFT) 
example show the behavior of spurious harmonics as a function of the signal-to - 
noise ratio (SNR). These results have strong implications for digital reception and 
signal analysis systems. At low SNRs , spurious harmonics decay exponentially on 
a log-log scale, and the resulting spectrum is white. As the SNR increases } the 
spurious harmonics figure prominently in the output spectrum. A useful expression 
is given that roughly bounds the magnitude of a spurious harmonic as a function 
of the SNR. 


I. Introduction 

This work was inspired by consideration of a 2-million 
channel spectrum analyzer built by the Digital Projects 
Group of the Communications Systems Research Sec- 
tion [1]. This spectrum analyzer is a prototype of a larger 
system that will be used in the sky-survey portion of 
the Search for Extraterrestrial Intelligence (SETI) project. 
After computer simulations were performed, 8-bit input 
quantization was observed to pose the greatest limitation 
to the dynamic range of the spectrum analyzer. This is 


1 The author is also a graduate student in Electrical Engineering at 
the California Institute of Technology. 


because quantization is a nonlinear process that generates 
spurious harmonics in the spectrum of the quantizer out- 
put. 

Previous work by Bennett [2] considered the spectra of 
quantized signals when the system input has “energy uni- 
formly distributed throughout a definite frequency band 
and with the phases of the components randomly dis- 
tributed.” Hurd [3] developed an expression for the cor- 
relation function of a quantized sine wave plus Gaussian 
noise and examined the case where the input noise spec- 
trum is rectangular narrow-band and the signal-to-noise 
ratio (SNR) is small. Quantization error spectra are most 
commonly assumed to be white [4]. This article derives an 
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expression for the spectrum of a quantized sine wave plus 
white Gaussian noise. An SNR transition region where the 
spectrum goes from being filled with spurious harmonics 
to white is presented. This transition is due to the dither- 
ing effect of the input Gaussian noise. A rule of thumb 
is given bounding the size of a spurious harmonic as a 
function of the SNR. Implications for digital reception and 
signal analysis systems are considered. 

II. Power Spectrum of the Quantization Error 

Consider the quantizer system in Fig. 1 with input x 
and output y. One can write: 

y - Q[x] = x - e (1) 

where Q[ ] is the quantization operator and e is the quan- 
tization error. When Q[ ] is a uniform mid-tread symmet- 
ric quantizer with a staircase input-output relation as in 
Fig, 2, e can be expressed as a sawtooth function of x as 
in Fig. 3. Assuming an infinite quantizer (or equivalently, 
no quantizer saturation), one can write a Fourier series 
expansion for e{x) as in [5]; 



fO 


Now consider the system in Fig. 4 where the input is 
A sin(u?o t + <f) + z{t ) and 2 r(t) is zero mean, Gaussian noise 
with variance er 2 . The continuous-time signal e(t) can be 
written as 


«(*) = 

Lib exp ^^^i(>lsin(wo< + ^) + «(0)) 

**o° 

(3) 

The autocorrelation function of e(t) is defined as 

R e (t,t + r) = E{e(t)e(t + r)} (4) 

Returning to Fig. 4, one sees that x[n] = e(nT) is a 
discrete-time random process. The autocorrelation func- 
tion of x[n] can be expressed as 

R x [n , n + k] = E{x[n]x[n -f k ]} 

= E{e(nT)e(nT + kT)} = R e (nT, nT + kT) 

(5) 

If the phase <f> of the input sinusoid is a random variable 
uniformly distributed between 0 and 27r, can be com- 

puted as 

— If 2 * 

= — / R x [n,n + k]d<f> (6) 

zir Jq 

Leaving the details to Appendix A, one obtains 


R x [k] 


(-1)' (1 , 

CO 

Y. S n exp(jmj 0 Tk); 


(¥) 


n — oo 
n odd 


k = 0 

k^O 


(7) 


where 


and J n (x) is the nth-order Bessel function. The power 
spectrum of the discrete-time random process x[n] is 



( 8 ) 


OO 

S(w) = £ ^lk]e juTk (9) 

k = - oo 
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Finally, the power spectrum can be expressed as 


oo 

S(u>) = N q + ^2 2irS n 6((ujT — nwoT) mod 2 tt) (10) 

n — oo 
n odd 


where 



n odd 


and S n is defined in Eq. (8) (details of this derivation 
are left to Appendix A). The Nq term in Eq. (10) rep- 
resents a frequency-independent quantization noise floor 
(see Appendix A). The infinite sum in Eq. (10) specifies 
the phase-averaged magnitude and location of all spurious 
harmonics in the frequency domain. Depending on the 
value of u>qT , spurious harmonics can be spread through- 
out the frequency domain or lie concentrated at only a few 
frequencies. This complicates digital reception of weak 
signals in the presence of stronger interferes. While a 
strong interferer may be easily identified and filtered out, 
the spurious harmonics generated by a strong interferer 
would require more complex filtering techniques. A spec- 
trum (taken from a discrete Fourier transform) with spurs 
is shown in Fig. 5. 

When a > A, 


mediately obvious how the spectrum will appear. For this 
reason, an example involving a discrete Fourier transform 
(DFT) is presented in the next section. 


III. DFT Example 

This example provides a quantitative analysis of the 
manner in which spurious harmonics are dithered due to 
additive white Gaussian noise. In particular, an SNR re- 
gion where the magnitude of a spurious harmonic decays 
exponentially on a log-log scale is presented. Consider the 
DFT of the signal x[n] in Fig. 4. One can write 


N-\ N-\ 


£{Pa*]| 2 } = iEE £Mn]*M}e _;2 T~ r) (11) 


n = 0 r = 0 


For the purpose of exposition, the u>qT product is chosen 
to equal 7 t/ 8, and the value of k is 3./V/16. Thus, a case is 
analyzed where the infinite number of spurious harmonics 
are aliased into only a few frequency bins. In this exam- 
ple, the phase is not treated as a random variable as was 
done in the previous section. Treating the phase as a con- 
stant allows a more general phase-dependent solution to 
be obtained. 


By observing bin 3 N/16, one is examining the spectral 
sample that contains the spurious harmonic specified by 
an arrow in Fig. 5. The spectrum in Fig. 5 was generated 
using the following parameters: A = 0.5, <f> = 0.8147576, 
N = 1024, A = 1/127, and SNR = 10 log 10 (A 2 /2<x 2 ) = 
50 dB. 


A 2 

= — + O 



— 2tt 2 <7 2 

A 2 


-|- < spur term > 


Leaving the details to Appendix B, one can evaluate 
Eq. (11): 


where the <spur term> is a delta function with weight zero y = £{|A[3A716]| 2 } = Nqi + N Q2 + N Q3 + Sq (12) 

or weight <9(exp(-47r 2 <r 2 /A 2 )), depending on whether or 

not a spur was located at that frequency. These results 

are consistent with those in [5]. where 


With the exception of the A 2 / 12 term, all the compo- 
nents in the power spectrum in Eq. (10) are dependent on 
the ratio <r/A. When this ratio is much greater than 1, the 
resulting spectrum is essentially A 2 /12 and white. In this 
manner, the input Gaussian noise has a dithering effect 
on the spectrum. For other values of <r/A, it is not im- and 


A 2 

Nqi = UN 


(13) 
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N, 


Q 2 


A 2 * “ < ! )' exp ( A 2 ,2 ) [ [2irAQ m l\ 

= Te^v ee — — p r — s — ; 


(2*AZ m l\] 

-f COS 7 • ! 

V A 




AT — ^ 

g3 “ 16n 2 N E 

m = 0 


(-l)'exp(-^ 2 ) 


sin 


/=i 


^ 2 t rAe m / j 


+ 


- (-D'exp(-^) 

S E E s i n 


E 

/-I 




2n 




oo (-l)'exp(-^E/ 2 ) 

E L^__^ F ( 16 . )3l 16,^) 

f=l v / 




J 


(14) 


with F( ), 0 m , and E m as defined in Appendix C. 

The quantities in Eq. (12) are easily computed. In fact, 
for the SNR s of interest, only a small number of terms is re- 
quired to adequately represent the infinite sums over / (see 
Appendix D for more details). The spectral sample Y con- 
sists of two types of components: a white-noise-floor term 
and a spurious-harmonic term. As seen in Appendix B, the 
values of Nq t , Nq 7 , and Nq s are independent of bin num- 
ber (frequency), and thus represent a white-noise floor. 
As a check, these quantities obey the 1/N (where N is the 
number of points in the DFT) processing rule for white 
noise in a DFT. The value Sq represents the spurious har- 
monic component of the spectral sample. This term is 
independent of the DFT size N and is responsible for the 
spectral sample becoming large when the SNR is high. In 
effect, this limits the dynamic range of the digital spec- 
trum analyzer. 

Figures 6(a) and (b) show how the value of the spectral 
sample changes as the SNR varies. The units are decibels 
relative to the carrier (dBc). The following parameters 
are held constant: A = 0.5, <j> — 0.8147576, N =■ 1024, 
and A = 1/127. The quantizer step size A was chosen 
to simulate an 8-bit input quantizer (A « 2“ B + 1 where 
5 = 8). An experimental curve (simulation) is presented 
with the theoretical curve to show the excellent agreement. 
For each SNR used in the experiments, 10,000 spectra were 
accumulated and rescaled. 


The horizontal line in Fig. 6(a) indicates the value of 
Nq x . This would be the value of the spectral sample un- 
der traditional quantization error assumptions [4]. Below 
40 dB SNR (for the 8-bit input quantizer), the value of 
the spectral sample Y reduces to essentially Nq 1 . Above 
70 dB SNR, the value of the spectral sample does not 
change by more than a few decibels. The transition region 
in Fig. 6(a) coincides with a j A approaching and exceeding 
unity. 

Figure 7 plots Sq , the spurious harmonic portion of the 
spectral sample, as a function of SNR. The units are deci- 
bels relative to the carrier (dBc) with the same conditions 
as above. When <t/A approaches unity, this value decays 
exponentially on the log-log scale. From Eq. (14), one can 
estimate the behavior of Sq when a /A exceeds unity to 
obtain a rule of thumb. Assume that the first term in 
the infinite sum over / in Eq. (14) dominates, and bound 
F ( ) by its maximum value (from Appendix C, \F( )| < 1). 
This provides a rough bound on the maximum value of Sq , 
In particular, compute SQdBc = lOlogioC^Sg/A 2 ). Recall 
that SNR = 10 log 10(A 2 /2 <t 2 ). After manipulating, 

. { 2A 2 \ 207T 2 A 2 -snr 

SqdBc = 101og 10 J “ T^Io A 1 10 10 ^ 

This expression is easily inverted to express the SNR as a 
function of Sqabc* As seen in Fig. 7, the rule of thumb 
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in Eq. (15) nicely describes the behavior of Sq even when 
< 7 / A is less than one. Asymptotically, Eq. (15) levels off 
at high SNRs to a constant value, which is consistent with 
empirical observations. This expression should prove use- 
ful to system designers concerned with dynamic-range lim- 
itations imposed by input quantization. 

IV, Conclusions 

Spurious harmonics pose complex filtering problems for 
digital reception systems. These harmonics also limit the 


dynamic range of digital spectrum analyzers. Expressions 
have been obtained describing the spectrum of quantiza- 
tion error when the input is a noisy sinusoid. An exam- 
ple involving a DFT has provided quantitative information 
about the behavior of spurious harmonics in the frequency 
domain as a function of the SNR. The input Gaussian noise 
dithers the output spectrum when the ratio <x/ A exceeds 
1, where a 2 is the noise variance, and A is the step size 
in the input quantizer. A useful rule of thumb has been 
derived that roughly bounds the magnitude of a spurious 
harmonic as a function of the SNR. 


Acknowledgment 


Simulations in this report were done in part on the JPL-CRAY. The author 
thanks the Supercomputer User Group for their assistance. 


159 


Fig. 1. Quantizer system diagram. 




Fig. 4. Generation of discrete-time quantization error. 



DISCRETE FOURIER TRANSFORM BIN 

Fig. 5. DFT spectrum, 50-dB SNR Input. 


Fig. 3. Quantization error input-output relation. 






Appendix A 

Derivation of the Power Spectrum 


From Eq. (4), first consider r = 0: 


A 2 


oo oo 


«.('.')= -j3 E E 


(_!)<+* 


4tt 2 ^ ^ Ik 

* = — OO /= — OO 
*5*0 


x exp 


/i27r>l 

V A 


(/ + fc) sin(w 0 < 4- <£)) 


x £ 


MiMt*)*,)} 


a2 00 1 a 2 00 

s? ,E ?+«? E c-»r 

^5*0 m ;*0 

x “ p (- 5 - S1 ” <"“* + *>J , E j(73S) 

?5*0 


Noting that the sum involving l// 2 evaluates to tt 2 /3 [7] 
and (after partial fraction manipulation) that the sum in- 
volving !//(/ — m) evaluates to 2/m 2 , one can write 


Using the characteristic function of a zero-mean Gaus- 
sian random variable z with variance <r 2 , it is known that 
i?{exp(jaz)} = exp(a 2 /2er 2 ) [6]. Therefore, 




A 2 

4tt 2 


00 oo 

E E 

* = — 00 /= — co 
*j *0 ^ 5*0 


Ik 


x exp 



(/ + £)sin 



x exp 


/ 2n' 1 <r 2 

V A 2 " 


n+k) 


A 2 

R'(t,t) = j 2 


A 2 ,22, (~l) m exp 8 in ( w 0 t + <£)) 

+ W E 


m= — 00 
m ;*0 




27T 2 <T 2 2^ 

Xexp V A*-™) C^' 1 ) 

Now consider r ^ 0: 

+ r) = 


Changing indices so that / + k = m and translating the 
condition k ^ 0 into m ^ l, one obtains 


a 2 “ ~ (_i y+t 

45T 2 E 2^ /jfc 

* = — 00 /= — oo 
*5*0 /5*0 


R e (t,t) = 

a 2 ~ ~ (-ir 

4a- 2 E 2^ 

m= — oo /= — oo v 7 

m ^0 * 5*0 


x exp f ■ [fc sin (u>ot -f <j>) *f / sin ( ujqt + <f>)] 

x E | exp (kz(t) + lz(t + r))^ j 


( j27rAm . , A f 2 tt 2 <7 2 2 \ 

x exp i — — — sin (wot + <j>) J exp I- nr 1 

Separating the case when m = 0 and rearranging the order 
of summation, 


Since the input noise is white, z(t) and z(t + r) are in- 
dependent random variables (r ^ 0), and the expectation 
of the product becomes the product of the expectations. 
When one uses the same characteristic function method 
detailed above, 
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A 2 ~ ~ (-1)'+* 

R e (t,t + r) — 4ff2 X/ ik 

t=-oo !=- oo 
ij50 

/ 2tt 2 (t 2 2 2 A 
x exp ^ &T~( l + m ) J 

x exp sin(u;o^ 4* <£) 

-f / sin(wo< 4 wot 4 <£)]^ (A-2) 

As indicated in Eq. (5), to obtain the discrete-time au- 
tocorrelation function R x [n } n4*]> replace t with nT and r 
with kT in Eqs. (A-l) and (A-2). The evaluation of Eq. (6) 
involves changing the order of integration and summation 
until only terms involving the phase <j> are inside the in- 
tegral. It is useful at this time to use the Jacobi-Anger 
formula [7]: 

oo 

e J ' I8in * = ^2 

p= — 00 


00 A 2 
R x [k] = ~ 2 (-l) p exp(-jpw 0 Tfc)^2 





\ 

X 

oo 

E 

tm. 130. exp ( 

f 2*V 2 2 \ 

K A2 ) 


/= — oo 
\ 1*0 


/ 


Consider the sum over L When p is even, this sum will 
be zero since J p (x) = (-l) p J p (-x). When p is odd, the 
sum over / is an even function of /, and one can reduce the 
double-sided infinite sum to a single-sided infinite sum. By 
changing the index from p to — n, one obtains the second 
half of Eq. (7). 

In evaluating Eq. (9), one can write 


S(w) = fl x [0] + ]T R x [k]e’ uTk 

Jb= — oo 
k*0 

Note that the first term, fl x [0], is independent of the fre- 
quency, w. Using the notation in Eq. (7), 


One can now evaluate 


TUo] = 


A 2 A 2 
l2 + 2r 2 



(- 1 )’ 

P 


* £ J - (rr) eW " 

\ ' p=— OO ^ ' 

*kj e “‘* a * 

0 

The integral above will be 1 when j> = 0 and 0 otherwise. 
Noting that the resulting expression inside the above sum 
is an even function of /, one obtains the first half of Eq. (7). 
In evaluating when k ^ 0, one again employs the 

Jacobi-Anger formula and interchanges the order of sum- 
mation and integration so that the appropriate terms from 
Eq. (A-2) are inside the integral. Using the orthogonality 
of exponentials and noting that J p (x) — (— 1 ) p J-p(x), one 
obtains 


S(u > ) =7k[0] 


4 


Y Y Sn exp(jnw 0 Tfc) exp(jwTi) 

kzz-oo 71-00 
i n odd 


oo 

- E s " 

n — oo 
n odd 


Note that the third term above is independent of fre- 
quency. Writing out the expression for Rg[ 0], one arrives 
at Eq. (10) after noting 


oo oo 

y] e* xh = 2ir ^ 6(x — 27rfc) = 2 ir6(x mod 2 tt) 

Jt=-oo fc=-00 


(A- 3) 


where S ( ) is the Dirac delta function. 
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Appendix B 

Derivation of the DFT Problem 


From Eq. (11), one can write 


y = E{\x[k]\ 2 } 


Nq 3 = 


f h: 

2n 2 N ^ m 2 

m = — oo 
m^O 


= 771 X) £{*[»]*[»]} 


n=r 0 


' ' r>— _ r*-\ V / 


p= — OO 


+ w £ E ®(«w»w)«p ( ,,2, y r) ) 

n — n n \ / 


n=0 r = 0 
r^n 


(B-l) 


Consider the first term above. Note that it is independent 
of the bin number k (i.e., independent of frequency). Re- 
calling the results of Eqs. (5) and (A-l), rearranging the 
order of summation, and using the Jacobi- Anger formula 
leaves one with 


The sum over p above is evaluated in Appendix C. Insert- 
ing the result of Appendix C, interchanging the order of 
summations, and using Euler’s rule yields the second part 
ofEq. (13). 

Now return to the second part of Eq. (B-l): 
w = E X) ex P ( ~ ;27r ^ n ~ r) ) 

n=0 r=0 ' ' 


N - 1 


) ! Re(nT, nT) — 


n= 0 


A 2 

12 N 


w T - «P (=^ 

n=0 r=0 ' 


- r ) 


r^n 


A 2 ^ 

*1* n *> xr 


2ir 2 N 


m = -oo 
m^O 


(- 1 )" 

m 2 


where R e ( ) is defined in Eq. (A-2). Evaluating this fur- 
ther, 


x -»(- ! sr"0 £ ^(t 2 )*** 

' ' n— — ro ' * 


P= — oo 


1 JV_1 
i V — ^ Jpm 

x — 2 ^ e * 


TV 2 


n=0 


Recall that u>qT = 7t/ 8 and = 3 jV/ 16 in this example. 
The sum over n is nonzero only when p = 0 mod 16. So, 


N— 1 


■^2 X] R'i n T,nT) = 


n=0 


where Nq , = A 2 /12iV, and 


W = 


A 2 


OO OO 


-4 D -^ ex p 


/=: — oo m = — oo 
1^0 m*0 


(-^ p 1 +■»■)) 


jj? S 51 “P [I rip (w 0 Tn + 


n=0 r=0 


4- m sin (o) 0 T r + <£)|j exp — — ) 

1 v-^ 1 f.2irA(l + m) N 

-Jpl^ ex P ( 3 a sm ( u '° Tn + *) 

n=0 ' ^ 
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Note that the second term inside the big brackets is in- 
dependent of the bin number k. Further separating the 
terms inside the big brackets, one can write 


W = Sq + Nq, (B-2) 


X ex p(~ 2 ^ + m2 )) 

[ JV-1AT-1 . 

■m H exp ( sin ( w ° rri + 

n= 0 r= 0 V 


Consider the Nq s term (this is the term independent of 
k). Employing the Jacobi- Anger formula and evaluating 
the sum over n as before, 


/ / — j2nk(n — r)\ 

+ msin(w 0 7Y 4* </>)]) exp ( — ) 


No. = 


A 2 

Att 2 N 


oo oo 


E E 

Izz — oo m — — oo 
1^0 0 


(_!)(»+"•) exp (- 1^1 (/2 + m 2 )) 

hn 


Call the sum in the big brackets V , and employ the Jacobi- 
Anger formula. Using the values of u^qT and k from Sec- 
tion III, and assuming that N is a power of 2 (N > 16), 
the value of V reduces to 


x E +"*))«'”** 

fc = - OO ^ 

The results of Appendix C can be used to evaluate the 
infinite sum over k. Again, after manipulating the defini- 
tion of F( ) in Appendix C as detailed above, one finally 
obtains the third part of Eq. (13). 


V = 


C 

x ( E J—.* (t 2 ) ^ h ) 

\q=- oo ' ' / 


Finally, consider Sq from Eq. (B-2): 


A 2 ^ A (-l) ,+m 

S Q--JZ2 2^ 2^ 


1 =i — oq m = — oo 
;^0 m? 0 


/m 


Since J n (x) = (-l) n J„(-x) and ±3 + 16p is odd, one 
can change the double-sided infinite sums to single-sided 
infinite sums. By changing the index on q to -q, recall- 
ing that J- n (x) = (-1)" J n (x), and using the formula in 
Appendix C, one gets Eq. (14). 
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Appendix C 

Evaluating an Infinite Bessel Sum 


Define 


oo 

F(x,y,z,t)= £ e"‘J, +,*(<) 

fc=-oo 


the remainder of this analysis assumes (for simplicity) that 
no delta functions lie on the limits of integration. For the 
DFT example in this article, this condition is satisfied. 
Using the sifting and scaling properties of delta functions, 
one finally obtains 


where x and y are integers and z / 0. Using the integral F(x } y> z } i) = 
definition of the Bessel function of integer order [7], 


F{.x,y,z,t) = 


U/2J-1 


^ exp (j (<0 m + y0 m )) + exp (j (t= m - y£ m ))j 


m=0 


oo -(-y+zk) } 

e ixkJ — /cos (2 ,6 + zkO)e> tQO *Wde 

*=-oo 0 

After changing the order of integration and summation, 
one can apply Euler’s rule to express cos( ) in terms of 
complex exponentials. Next, equate the infinite sums of 
complex exponentials to infinite sums of Dirac delta func- 
tions as in Eq. (A-3). Only a finite number of delta func- 
tions will remain inside the limits of integration. Note that 


where 



and 0m — cos(0 m ) and — cos(4m). 
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Appendix D 

Truncating an Infinite Sum 


This appendix considers bounds on the error introduced 
by truncating the following infinite sum: 


oo -a 3 / 3 


= E : 


(=1 


l n 


■m 





■m+ E 

J=L+ 1 



m 


=z x + error 

where n > 1 and |F(/)| < 1. The error can be bounded as 
follows: 

oo -a 3 ! 3 °° 

|error| < ^ , n |F(/)[ < E ex P( _Q2 * 2 ) 

I=i,+ 1 1=L+ 1 


Let a 2 = 1/2 tr 2 . Then, 


OO 

|error| < y/2na 2 y] 
f=i+i 


V2 


Z 6 




The above sum can be visualized as the integral from 
x — I to a: = oo of a discontinuous step-like function 
whose value is constant over the interval between two ad- 
jacent integers. The value over an interval is equal to the 
Gaussian density evaluated at the right-most portion of 
the interval A little thought will show that this integral 
is strictly less than the integral of a Gaussian distribution 
from x = L to x = oo. Therefore, 


|error| < ^^erfc(al) 
2a 


where erfc( ) is the complementary function defined in [7], 
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